Демонстрация рядов

ts_data3 <- ts(read.csv("ts3.txt", header = TRUE, as.is = FALSE, sep = ','))
plot(ts_data3) 

ts_data7 <- ts(read.csv("ts7.txt", header = TRUE, as.is = FALSE, sep = ','))
plot(ts_data7)

Ряд 3 не колеблется вокруг 0, значит в модели есть константа, ряд нестационарный.

Ряд 7 колеблется вокруг 0, значит в модели нет константы, предположим, что ряд стационарный.

Ряд 3

Подбор модели

acf(ts_data3)

pacf(ts_data3)

adf.test(ts_data3)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 0.7251   0.853
## [2,]   1 0.1004   0.673
## [3,]   2 0.0881   0.669
## [4,]   3 0.0717   0.665
## [5,]   4 0.0314   0.653
## [6,]   5 0.0584   0.661
## [7,]   6 0.0642   0.662
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -2.70  0.0788
## [2,]   1 -2.07  0.3017
## [3,]   2 -2.05  0.3096
## [4,]   3 -2.11  0.2826
## [5,]   4 -2.06  0.3034
## [6,]   5 -2.19  0.2530
## [7,]   6 -2.20  0.2492
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -3.20  0.0876
## [2,]   1 -2.29  0.4550
## [3,]   2 -2.27  0.4614
## [4,]   3 -2.34  0.4338
## [5,]   4 -2.28  0.4593
## [6,]   5 -2.42  0.4007
## [7,]   6 -2.44  0.3923
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Смотрим на p-value для модели с константой (2 и 3), p-value во всех случаях больше 0.05, считаем, что присутствует единичный корень, а значит продифференцируем ряд.

ts_data3_diff <- diff(ts_data3, differences = 1)
plot(ts_data3_diff)

adf.test(ts_data3_diff)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 -12.73    0.01
## [2,]   1 -11.70    0.01
## [3,]   2 -10.39    0.01
## [4,]   3  -9.51    0.01
## [5,]   4  -8.98    0.01
## [6,]   5  -8.69    0.01
## [7,]   6  -8.07    0.01
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -12.73    0.01
## [2,]   1 -11.70    0.01
## [3,]   2 -10.39    0.01
## [4,]   3  -9.51    0.01
## [5,]   4  -8.99    0.01
## [6,]   5  -8.69    0.01
## [7,]   6  -8.07    0.01
## Type 3: with drift and trend 
##      lag    ADF p.value
## [1,]   0 -12.98    0.01
## [2,]   1 -11.96    0.01
## [3,]   2 -10.65    0.01
## [4,]   3  -9.75    0.01
## [5,]   4  -9.26    0.01
## [6,]   5  -8.98    0.01
## [7,]   6  -8.36    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Теперь ряд колеблется вокруг 0. И при этом p-value в модели без сдвига случаях меньше 0.05, считаем, что нет единичного корня и \(d\) = 1, так как хватило одного дифференцирования.

acf(ts_data3_diff)

pacf(ts_data3_diff)

На графике acf видим убывание, похожее на \(e^x\) при этом на графике pacf значима только первая корреляция, считаем \(p\) = 1 и \(q\) = 0. (Если было бы \(q > 0\), то acf быстрее бы убывало, при этом \(p\) определяем как количество значимых корреляций на pacf)

Получили модель ARIMA(1,1,0). Положительное \(\varphi \approx 0.7\).

auto.arima

aa3 <- auto.arima(ts_data3, test = 'adf')
aa3
## Series: ts_data3 
## ARIMA(2,1,1) 
## 
## Coefficients:
##           ar1     ar2     ma1
##       -0.2237  0.6776  0.9585
## s.e.   0.0713  0.0578  0.0610
## 
## sigma^2 = 1.068:  log likelihood = -1450.93
## AIC=2909.86   AICc=2909.9   BIC=2929.49

Выходит совсем другая модель :)

Подбор параметров

model3_my <- Arima(ts_data3, order = c(1, 1, 0))
checkresiduals(model3_my)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 8.2748, df = 9, p-value = 0.5067
## 
## Model df: 1.   Total lags used: 10
checkresiduals(aa3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)
## Q* = 6.574, df = 7, p-value = 0.4745
## 
## Model df: 3.   Total lags used: 10

В обеих случаях остаток- белый шум \(\sim N(0,\sigma^2)\).

Прогноз

Построим прогноз для обеих моделей (но мне больше нравится моя модель, так как в ней меньше параметров).

Отрежем от ряда 65 точек (после этой точки резкое убывание) и построим прогнозы по моделям, подогнанным под обрезанный ряд.

ts_data3_fc <- ts_data3[1:935]
plot(forecast::forecast(Arima(ts_data3_fc, order = c(1, 1, 0)), h = 65))
lines(936:1000, ts_data3[936:1000], col = 'red')

plot(forecast::forecast(Arima(ts_data3_fc, order = c(2, 1, 1)), h = 65))
lines(936:1000, ts_data3[936:1000], col = 'red')

Прогнозы очень близки и при этом из-за резкого убывания графика рельное продолжение попадает только в 95% интервал.

Ряд 7

Выбор модели

adf.test(ts_data7)
## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -7.50    0.01
## [2,]   1 -8.76    0.01
## [3,]   2 -8.12    0.01
## [4,]   3 -7.59    0.01
## [5,]   4 -7.72    0.01
## [6,]   5 -7.49    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -7.58    0.01
## [2,]   1 -8.88    0.01
## [3,]   2 -8.24    0.01
## [4,]   3 -7.72    0.01
## [5,]   4 -7.88    0.01
## [6,]   5 -7.68    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -7.60    0.01
## [2,]   1 -8.91    0.01
## [3,]   2 -8.27    0.01
## [4,]   3 -7.75    0.01
## [5,]   4 -7.91    0.01
## [6,]   5 -7.70    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

p-value для модели без константы во всех случаях меньше 0.05, не принимаем гипотезу о наличии единичного корня (процесс стационарный). В модели считаем \(d\) = 0.

acf(ts_data7)

pacf(ts_data7)

На графике acf видим убывание, похожее на \(e^x\cdot\cos{x}\) и на графике pacf значима две корреляции, считаем \(p\) = 2 и \(q\) = 0.

Получили модель ARIMA(2,0,0).

auto.arima

aa7 <- auto.arima(ts_data7, test = 'adf')
aa7
## Series: ts_data7 
## ARIMA(2,0,0) with zero mean 
## 
## Coefficients:
##          ar1      ar2
##       0.9625  -0.2080
## s.e.  0.0437   0.0438
## 
## sigma^2 = 0.9723:  log likelihood = -701.99
## AIC=1409.97   AICc=1410.02   BIC=1422.61
checkresiduals(aa7)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0) with zero mean
## Q* = 3.48, df = 8, p-value = 0.9007
## 
## Model df: 2.   Total lags used: 10

Вау, совпало.

Прогноз

Отрежем от ряда 31 точку (после этой точки резкое возрастание) и построим прогноз по модели, подогнанной под обрезанный ряд.

ts_data7_fc <- ts_data7[1:469]
plot(forecast::forecast(Arima(ts_data7_fc, order = c(2, 0, 0)), h = 31))
lines(470:500, ts_data7[470:500], col = 'red')

Реальное продолжение ряда помещается в 95% интервал и колеблется вокруг спрогнозированной линии.

Сравнение прогнозов SSA, seasonal ARIMA и ETS

ts_data <- read.csv("MRTSSM4453USN.csv", header = TRUE, as.is = FALSE, sep = ',')[,1:2]
ts_data[,1] <- c(0:359)
ts1 <- ts_data[1:336,]
ts1 <- ts(ts1[,2], start = c(1992, 1), frequency = 12)
ts1s <- ts(ts1[1:324], start = c(1992, 1), frequency = 12)
ts1s2 <- ts(ts1[1:312], start = c(1992, 1), frequency = 12)
ssa.result_1 <- ssa(ts1s, L=324/2)
ssa.result_2 <- ssa(ts1s2, L=312/2)

Для одного периода:

auto.arima(ts1s)
## Series: ts1s 
## ARIMA(3,1,2)(0,1,2)[12] 
## 
## Coefficients:
##           ar1     ar2     ar3      ma1      ma2     sma1     sma2
##       -0.3078  0.1964  0.3714  -0.5782  -0.3112  -0.2548  -0.1585
## s.e.   0.1799  0.0944  0.0719   0.1854   0.1438   0.0608   0.0530
## 
## sigma^2 = 5872:  log likelihood = -1788.86
## AIC=3593.73   AICc=3594.2   BIC=3623.65
bf_v <- forecast::forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')

bvec_1 <- rmse(ts1[325:336], bf_v$mean)
ar_1 <- forecast::forecast(Arima(ts1s, order = c(3, 1, 2), seasonal = c(0, 1, 2)), h = 12)
plot(ar_1)
lines(ts1, col = 'black', t='l')

fora_1 <- rmse(ts1[325:336], ar_1$mean)
ets_1 <- forecast::forecast(ets(ts1s), h = 12)
forets_1 <- rmse(ts1[325:336], ets_1$mean)
plot(ets_1)
lines(ts1, col = 'black', t='l')

bvec_1
## [1] 85.11311
fora_1
## [1] 97.87564
forets_1
## [1] 120.1644

Лучший результат– SSA с MSE = 85.11.

auto.arima(ts1s2)
## Series: ts1s2 
## ARIMA(3,1,2)(0,1,2)[12] 
## 
## Coefficients:
##           ar1     ar2     ar3      ma1      ma2     sma1     sma2
##       -0.2954  0.2044  0.3700  -0.5788  -0.3164  -0.2615  -0.1522
## s.e.   0.1867  0.0972  0.0742   0.1921   0.1480   0.0632   0.0546
## 
## sigma^2 = 5792:  log likelihood = -1717.72
## AIC=3451.44   AICc=3451.93   BIC=3481.04
bf_r <- forecast::forecast(ssa.result_2, groups = list(1:12), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')

brec_2 <- rmse(ts1[313:336], bf_r$mean)
ar_2 <- forecast::forecast(Arima(ts1s2, order = c(3, 1, 2), seasonal = c(0, 1, 2)), h = 24)
plot(ar_2)
lines(ts1, col = 'black', t='l')

fora_2 <- rmse(ts1[313:336], ar_2$mean)
ets_2 <- forecast::forecast(ets(ts1s2), h = 24)
forets_2 <- rmse(ts1[313:336], ets_2$mean)
plot(ets_2)
lines(ts1, col = 'black', t='l')

brec_2
## [1] 134.0593
fora_2
## [1] 107.2831
forets_2
## [1] 129.947

Лучший результат– seasonal ARIMA с MSE = 107.28.

MSSA

Построим два модельных ряда, первый– квадратичный тренд + синус с большим периодом + \(e^x\cdot cos(2\pi\omega n)\) + достаточно сильный шум; второй– квадратичный тренд + \(e^x\cdot cos(2\pi\omega n)\) с меньшей амплитудой + шум слабее.

a <- -1/600
b <- 12
c <- 5
n <- c(1:1008)
ts_mod1 <- ts(c*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100 + sin(2*pi*n/500)) + rnorm(1008, 0, 0.45)
ts_mod2 <- ts(2*c/3*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100) + rnorm(1008, 0, 0.15)
plot(ts_mod1, type = 'l')
plot(ts_mod1, type = 'l')
lines(ts_mod2, type = 'l', col = 'red')

Сравним результаты ssa и mssa:

ssa_signal_1 <- ssa(ts_mod1, svd.method = 'svd')
plot(ssa_signal_1, type = 'vectors')

ssa_signal_2 <- ssa(ts_mod2, svd.method = 'svd')
plot(ssa_signal_2, type = 'vectors')

ssa_signal_m <- ssa(cbind(ts_mod1, ts_mod2), svd.method = 'svd', kind = 'mssa')
plot(ssa_signal_m, type = 'vectors')

r1 <- reconstruct(ssa_signal_1,
                       groups = list(Trend = c(1, 4:6),
                                  Seasonality = c(2:3)))
r2 <- reconstruct(ssa_signal_m,
                       groups = list(Trend = c(1, 4:6),
                                  Seasonality = c(2:3)))
plot(r2, add.residuals = FALSE,
     plot.method = "xyplot",
     slice = list(component = 1), 
     screens = list('ts_mod1', 'ts_mod2'),
     col = 
       c("blue", "green"),
     lty = rep(c(1, 2), each = 2),
     scales = list(y = list(draw = FALSE)),
     layout = c(1, 2))

plot(r2, plot.method = "xyplot", add.original = FALSE,
     add.residuals = FALSE, slice = list(component = 2),
     col = 
       c("blue", "green"),
     scales = list(y = list(draw = FALSE)),
     layout = c(1, 2))

plot(r1$Trend + r1$Seasonality, t = 'l')

plot(r2$Trend[1:1008] + r2$Seasonality[1:1008], t = 'l')

plot(ts_mod1, t = 'l')
lines(r1$Trend, t = 'l', col = 'blue')
lines(r2$Trend[1:1008], t = 'l', col = 'red')

Восстановленные ряды не отличаются, тренды отличны только в начале при этом очень слабо.

На данном примере результаты очень близки, нельзя сказать, что какой-то из методов показывает результаты лучше.

2D-SSA

Взяли изображение с шумом.

img <- readJPEG("12.jpg", native =TRUE)
plot(0:2, 0:2, type='n')
rasterImage(img, 0, 0, 2, 2)

2-D SSA с окном (10, 10):

res_ssa <- ssa(img, kind = "2d-ssa", svd.method = "eigen", L = c(10,10))
plot(res_ssa, type = "vectors", idx = 1:20,
     cuts = 255, layout = c(10, 2),
     plot.contrib = FALSE)

plot(wcor(res_ssa, groups = 1:30),
     scales = list(at = c(10, 20, 30)))

r <- reconstruct(res_ssa, groups = list(Signal = c(1:4, 16)))
plot(r, cuts = 255, layout = c(3, 1))

В результате шум удален, возможно не полностью, попробуем другую длину окна (26,26):

res_ssa <- ssa(img, kind = "2d-ssa", svd.method = "eigen", L = c(26,26))
plot(res_ssa, type = "vectors", idx = 1:20,
     cuts = 255, layout = c(10, 2),
     plot.contrib = FALSE)

plot(wcor(res_ssa, groups = 1:30),
     scales = list(at = c(10, 20, 30)))

r <- reconstruct(res_ssa, groups = list(Signal = c(1:10)))
plot(r, cuts = 255, layout = c(3, 1))

Шум удаляется, контуры чуть более четкие, но понадобилось большее число компонент.